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ABSTRACT 

We study the radiation-driven warping of accretion discs in the context of X-ray 
binaries. The latest evolutionary equations are adopted, which extend the classical 
alpha theory to time-dependent thin discs with non-linear warps. We also develop 
accurate, analytical expressions for the tidal torque and the radiation torque, including 
self-shadowing. 

We investigate the possible non-linear dynamics of the system within the frame- 
work of bifurcation theory. First, we re-examine the stability of an initially flat disc to 
the Pringle instability. Then we compute directly the branches of non-linear solutions 
representing steadily precessing discs. Finally, we determine the stability of the non- 
linear solutions. Each problem involves only ordinary differential equations, allowing 
a rapid, accurate and well resolved solution. 

We find that radiation-driven warping is probably not a common occurrence in 
low-mass X-ray binaries. We also find that stable, steadily precessing discs exist for 
a narrow range of parameters close to the stability limit. This could explain why so 
few systems show clear, repeatable 'super-orbital' variations. The best examples of 
such systems. Her X-1, SS 433 and LMC X-4, all lie close to the stability limit for a 
reasonable choice of parameters. Systems far from the stability limit, including Cyg 
X-2, Cen X-3 and SMC X-1, probably experience quasi-periodic or chaotic variability 
as first noticed by Wijers & Pringle. We show that radiation-driven warping provides a 
coherent and persuasive framework but that it does not provide a generic explanation 
for the long-term variabilities in all X-ray binaries. 
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1 INTRODUCTION 



light-curves (e.g. Gerend & Boynton 1976; Leibowitz 1984; 
Heemskerk & van Paradijs 1989; Still at al. 1997). 



A growing number of X-ray binaries are known to exhibit 
more or less periodic variability on time-scales much longer 
than their orbital periods (see Priedhorsky & Holt 1987; 
Smale & Lochner 1992; White, Nagase & Parmar 1995; Cor- 
bet & Peele 1997). The clearest examples of such 'super- 
orbital' periods can be found in Her X-1, SS 433 and LMC 
X-4. They have been known for a long time and are well 
documented. Different explanations have been proposed for 
the long periods of these systems but there seems to be a 
converging opinion, from both the observational and theo- 
retical viewpoints, that these periods represent the preces- 
sion of a tilted accretion disc (Katz 1973; Petterson 1975; 
Schwarzenberg-Czerny 1992; Wijers & Pringle 1999 and ref- 
erences therein). As well as periodically obscuring the X-ray 
source, a tilted disc casts an asymmetric shadow on the ir- 
radiated companion star which well reproduces the optical 



It has long been recognized that the precession of a 
tilted disc could be caused by the tidal force of the compan- 
ion star (Katz 1973). Consider the disc to be composed of 
a sequence of concentric circular rings, tilted with respect 
to the binary orbital plane. The rings behave gyroscopically 
owing to their rapid rotation. If they did not interact with 
each other, the tidal torque would cause each ring to pre- 
cess retrogradely, about an axis perpendicular to the binary 
plane, at a rate that depends on the radius of the ring, re- 
sulting in a rapid twisting of the disc. However, a fluid disc 
tends to resist this by establishing an internal torque be- 
tween neighbouring rings. This can be arranged so that the 
net torque on each ring is such as to produce a single, uni- 
form precession rate, allowing the disc to precess coherently. 
However, the internal torque can only be established if the 
disc becomes warped. This is accompanied by dissipation of 
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energy and the disc settles into the binary plane (Lubow & 
Ogilvie 2000; Bate et al. 2000). 

A mechanism is therefore required to excite the tilt con- 
tinuously. Radiation-driven warping has been identified as 
a plausible way to do this. If the outer part of the disc is 
warped, it intercepts radiation from the central point source. 
If this radiation is absorbed and re-emitted parallel to the lo- 
cal normal to the disc surface, the disc experiences a torque 
from the radiation pressure (Petterson 1977; Pringle 1996). 
The albedo of the disc is probably irrelevant since scattered 
photons leaving the disc will also have a distribution of mo- 
menta peaked around the local normal (see Appendix A4). 
A small tilt grows exponentially if it is twisted such that the 
line of nodes forms a leading spiral, and if the luminosity 
is sufficient to overcome dissipative processes that tend to 
flatten the disc (Pringle 1996). 

Such a mechanism would apply very well to X-ray bina- 
ries where the optical emission of the accretion disc is dom- 
inated by X-ray reprocessing (van Paradijs & McClintock 
1994). This does not preclude other driving mechanisms for 
warps, in particular wind torques (Schandl & Meyer 1994), 
but radiation-driven warping does have the advantage of re- 
lying on well-defined physics. 

The first question to be addressed is whether an ini- 
tially flat disc is unstable to radiation-driven warping. Pre- 
vious studies (Pringle 1996; Maloney, Begelman & Pringle 
1996; Maloney & Begelman 1997; Maloney, Begelman & 
Nowak 1998) concluded that discs in X-ray binaries are in- 
deed likely to be unstable. Although instructive, these stud- 
ies neglected self-shadowing of the disc and were not satis- 
fying in their treatment of the boundary conditions, which 
can be crucial to this problem. Furthermore, there was some 
uncertainty as to what equations to use. Recent progress 
(Ogilvie 1999, 2000) has led to a practical set of equations 
describing non-linear warps in thin discs. Although only one- 
dimensional, the equations are derived directly from the ba- 
sic fluid-dynamical equations in three dimensions, within the 
context of the alpha theory. 

The second question concerns the non-linear develop- 
ment of the instability. In particular, does it lead to steadily 
precessing discs as suggested by observations of some sys- 
tems? Recently, Wijers & Pringle (1999) investigated numer- 
ically the non-linear evolution of the instability and found 
a complex variety of behaviour. Their study was based on 
a direct integration of the earlier evolutionary equations 
of Pringle (1992). The calculations were restricted to low 
resolution owing mainly to the laborious treatment of self- 
shadowing. 

In this paper, we present a complementary study to 
that of Wijers & Pringle (1999). In Section 2 we present our 
formulation of the dynamical problem, based on the new 
evolutionary equations of Ogilvie (2000) . We also derive ac- 
curate, analytical expressions for the tidal torque and the 
radiation torque, including self-shadowing (see the Appen- 
dices). In Section 3 we re-examine the stability of an ini- 
tially flat disc to radiation-driven warping. We then search 
directly for non-linear solutions representing steadily pre- 
cessing discs and determine their stability (Section 4). All 
this can be achieved with high precision and resolution be- 
cause it involves solving only ordinary differential equations 
(ODEs). These solutions would be directly applicable to sys- 
tems such as Her X-1. Finally, in Section 5 we compare our 



results with observations. For the most part, we investigate 
generic properties in this paper, and leave the detailed mod- 
elling of individual systems to future work. 



2 DYNAMICS OF WARPED DISCS 
2.1 Background 

In general, the dynamics of a warped accretion disc in an 
X-ray binary is an extremely difficult problem of radiation 
magnetohydrodynamics in three dimensions. Furthermore, 
the evolution must be followed for many times longer than 
the characteristic orbital time-scale of the disc if one is to 
account for the observed behaviour. These considerations 
mean that an ab initio direct numerical treatment of the 
problem is completely unfeasible at present. 

Considerable progress has been made in recent years to- 
wards developing a reduced description of warped accretion 
discs that is only one-dimensional and follows the evolution 
on the long, viscous time-scale. This approach was intro- 
duced by Papaloizou & Pringle (1983) and developed by 
Pringle (1992). It is based on writing down plausible con- 
servation equations in one dimension for a disc composed 
of a continuum of thin circular rings that are in centrifugal 
balance but have varying inclinations. The rings interact by 
means of viscous torques. 

Although this approach is intentionally simplistic in ori- 
gin (Papaloizou & Pringle 1983), it has been shown that the 
same form of equations, with minor modifications, can be 
derived formally from the Navier-Stokes equation in three 
dimensions (Ogilvie 1999). This derivation also provides the 
non-trivial relation between the effective viscosity coeffi- 
cients appearing in the reduced equations and the small- 
scale effective viscosity that acts on motions with a length- 
scale comparable to the disc thickness. The resulting theory 
can then be seen as a non-linear extension of the existing 
linear theory of bending waves in viscous Keplerian discs 
(Papaloizou & Pringle 1983). 

To make this derivation possible, of course, certain as- 
sumptions must be made. The first is that the disc is thin, 
which is expected to be satisfied under a wide range of con- 
ditions. Indeed, the resulting theory is formally an asymp- 
totic theory for thin discs. The second is that the agent re- 
sponsible for angular momentum transport in accretion discs 
can be treated as an isotropic effective viscosity. Although 
it is now widely believed that this agent is magnetohydro- 
dynamic turbulence resulting from the magnetorotational 
instability (Balbus & Hawley 1998), our understanding of 
how the turbulence responds to imposed motions other than 
a Keplerian shear fiow is extremely limited (Abramowicz, 
Brandenburg & Lasota 1996). A first attempt to measure 
this response for motions of the type induced in a warped 
disc (Torkelsson et al. 2000) suggests that the assumption 
of an isotropic effective viscosity may be reasonable as a 
first approximation, and, in any case, no better model has 
been proposed. The third assumption is that the resonant 
bending-wave propagation found in inviscid, Keplerian discs 
may be neglected. This condition is satisfied if the viscos- 
ity parameter exceeds the angular semi-thickness of the disc 
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(a <; Ty/r).^] Further assumptions are of a more technical 
nature, but include the assumption that the fluid is poly- 
tropic. 

In a recent paper (Ogilvie 2000) this method was ex- 
tended to include non-trivial thermodynamics. Instead of 
assuming the fluid to be polytropic, one now accounts fully 
for viscous dissipation and radiative energy transport within 
the Rosseland approximation. The fluid is treated as an ideal 
gas, and a power-law opacity function, such as Thomson or 
Kramers opacity, is assumed. The resulting theory forms a 
closed system and may be considered as the logical exten- 
sion of the alpha-disc theory (Shakura & Sunyaev 1973) to 
time-dependent, warped accretion discs. In the absence of a 
warp, it reduces exactly to the usual alpha theory, except 
that the vertical structure of the disc is computed accu- 
rately by solving the differential equations, rather than by 
the usual algebraic approximations. 

2.2 Basic equations 

We consider a binary system with a compact star (or black 
hole) of mass Mi, about which the disc orbits, and a com- 
panion star of mass AI2, in a circular orbit of radius rb and 
angular frequency 

1/2 

f^b= rArii_^:ii^ . (1) 



+ Trad + Ttid 



(4) 



We adopt Cartesian coordinates {x, y, z) in a non-rotating 
frame of reference centred on the compact star, such that 
the orbit of the companion lies in the xy-plane and has a 
positive sense. 

The dynamical system is described by partial differen- 
tial equations in coordinates (r, i), where r is the distance 
from the central object and t is the time. At radius r, the 
disc rotates in centrifugal balance, with orbital angular ve- 
locity 0,{r) and specific angular momentum h{r) = r^Q.. The 
state of the disc is characterized by a surface density E(r, t) 
and a tilt vector l{r, t), which is a unit vector parallel to the 
local orbital angular momentum of the disc. The orbital an- 
gular momentum density is then L(r, i) — S/i I. Throughout 
this paper we take the angular velocity to be Keplerian, i.e. 
n = {GMi/r-y^^. 

The general forms of the conservation equations for 
mass and angular momentum are 



aE 1 a , ^, ^ 

and 



dL 1 a , ,, 

— + - — {rvL) 
ot r or 



IdG 
r dr 



(2) 



(3) 



respectively. Here v{r, t) is the mean radial velocity and 
27rG(r, i) is the internal torque in the disc. The terms S{r, t) 
and T(r, t) represent the sources (per unit area) of mass and 
angular momentum, respectively. If necessary, these terms 
should be considered as the mean sources averaged over the 
orbital time-scale. 

The total external torque consists of three parts, 

* Even if o ;$ H/r, resonant wave propagation is likely to be com- 
promised by non-linear effects, including a parametric instability 
(Gammie, Goodman & Ogilvie 2000). 



The first represents the angular momentum added to the 
disc along with the accretion stream. The radiation torque 
Trad, including the effects of self-shadowi ng, i s evaluated in 
Appendix A, where we find (cf. equation AlE) 



U 



127rrc 



■I X a. 



(5) 



Here L* is the luminosity of the central source (assumed 
isotropic) , c the speed of light and a a dimensionless vector 
depending on the shape of the disc. The tidal torque Ttidc, 
averaged over the binary or bit, is evaluated in Appendix B, 
where we find (cf. equation 1 



(6) 



3GM2Er% ,,, ,, ,^ , , 
Ttidc = —3 (e^ • 0(ez X Z) (IH ) . 

The final factor (not shown explicitly here) represents an 
advance over previous work; to omit it typically results in an 
error of 15% at r/vh — 0.3. Any other potential torques, such 
as self-gravitational, magnetic or gravitomagnetic torques, 
are neglected. Note that the acceleration of the origin of the 
coordinate system does not contribute a torque (Wijers & 
Pringle 1999). 



2.2.1 Internal torque 

Detailed analysis of the internal hydrodynamics of a warped 
disc (Ogilvie 1999, 2000) results in an expression for the 
internal torque of the form 

dl „ . dl- 



G = Ir^n^ [Qi I + Q2r-^ +Q3rlx—j, 



where 



pz dz 



(7) 



(8) 



is the second vertical moment of the density, and Qi , Q2 and 
Q3 are dimensionless coefficients. (In this definition only, z 
denotes the distance from the mid-plane of the warped disc, 
and the angle brackets denote azimuthal averaging.) The 
generalized alpha theory (Ogilvie 2000) provides a relation 
between X and E of the form 



X=Q5Cia^'"E 



l/7y.l0/7o-12/7 



n- 



(9) 



where Q5 is a further dimensionless coefficient, equal to 
unity for a flat disc, Ci is a constant, and a is the di- 
mensionless viscosity parameter.FI It is assumed here that 
the disc is Keplerian and gas-pressure-dominated, and that 
the opacity law is of the Kramers form k cx: pT~^'^. The 
value of the constant Cx is approximately 4.4 x 10^^ in CGS 
units (cm2°/''g-='/''s"^^/^) (Ogilvie 2000). Ahhough in the 
inner parts of the disc Thomson opacity may in fact dom- 
inate, and radiation pressure may be important, it will be 
seen that these inner parts play only a passive role in the 
dynamics, and equation (H) is sufficient for our purposes. 

T This equation is equivalent to the usual relation between the 
integrated viscous stress 'vS' and the surface density S in a flat 
disc. The coefficient Q5 reflects the thickening of the disc in re- 
sponse to increased viscous dissipation associated with motions 
induced by the warp. 
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The dimensionless coefficients Qi depend, most impor- 
tantly, on the shear viscosity parameter, a, and on the di 
mensionless amplitude of the warp, 

dl 



Mh 

= d(r - 

2-Kr ^ 



(14) 



1^1 = 



9 In] 



(10) 



There are much less significant dependences on the adiabatic 
exponent 7 and the bulk viscosity parameter ab- We assume 
that 7, a and Qb are constants, and set 7 = 5/3 and cvb = 
throughout. 

In this formulation, the internal torque is derived con- 
sistently from the basic fluid-dynamical equations in three 
dimensions, within the context of the alpha theory (i.e. the 
assumption that the effective dynamic viscosity is propor- 
tional to the local pressure). The torque responds dynami- 
cally to changes in E and increases in response to the en- 
hanced dissipation that results from a strong warp. This is 
different from the formulation of Wijers & Pringle (1999), 
who assumed the kinematic viscosity ;^ to be a fixed function 
of radius, independent of E and \'>p\. 

For small-amplitude warps, the coefficients have expan- 



Q. = Q,o + 0{\i!\^), 



(11) 



where 



Q30 = 



3a 
3(1 - 2a^ 



a(4 + a^ 



= 1. 



(12) 



2(4 + a2) ' 

We will consider first the case of linear hydrodynamics, in 
which the Qi are simply given the values Qio, independent 
of \ip\. Note that, even under this assumption, equations (GI) 
and (pi) do not form a linear system. We will also consider 
non-lmear hydrodynamics, in which the dependences of the 
Qi on \ip\ (and on 7 and Qb) are fully taken into account by 
solving a certain system of ODEs (Ogilvie 2000). Since only 
\ip\ varies within a given disc model, one can tabulate the 
functions (3i(|i/'|) beforehand, and simply interpolate during 
the course of the calculation. In our calculations the ampli- 
tude never exceeded \^\ — 1. 

2.2.2 Mass input 

In the case of a disc in a low-mass X-ray binary, the mass 
source 5* is associated with the accretion stream that origi- 
nates at the Li point on the surface of the Roche-lobe-filling 
companion star. When the disc is warped, the stream will 
intersect it at a point that depends on the instantaneous 
shape of the disc and on the binary phase (Schandl 1996) . 
In general, the stream reaches in to the vicinity of its cir- 
cularization radius re, which is much smaller than the outer 
radius ro of the disc. Although it is feasible to model this 
interaction to some extent within the framework of equa- 
tions (pl) and (bl), for simplicity, we follow Wijers & Pringle 
(1999) in assuming that the mass input is localized at a sin- 
gle radius which may be identified approximately as re. The 
added mass has the same specific angular momentum as a 
circular Keplerian orbit of radius re and aligned with the 
binary plane. We then have 



S = 



oir - 

27rr ^ 



(13) 



where M is the steady rate of mass supply. We defer a more 
accurate treatment to future work. 



2.2.3 Boundary conditions 

At the inner radius r = ri of the disc, there are several dif- 
ferent possible situations, depending on the nature of the 
central object. In the case of a black hole, the disc is likely 
to end (effectively) near the marginally stable circular orbit, 
where it makes a transition to a rapidly plunging spiral flow. 
At such a point, the surface density and internal torque fall 
essentially to zero (although see Krolik 1999 for an alterna- 
tive viewpoint). This possibility may also occur in the case 
of a weakly magnetized neutron star; but if not, the disc is 
likely to extend up to a boundary layer on the stellar sur- 
face. Finally, in the case of a strongly magnetized neutron 
star, the disc is likely to be truncated at a magnetospheric 
boundary. In this paper we apply the simple boundary con- 
ditions G = and dl/dr = at the inner radius. (These are 
independent, since the first already implies T — 0.) These 
conditions state that the internal torque vanishes at the in- 
ner radius, and that the disc is locally flat; however, the 
inclination of the inner edge is allowed to vary freely with- 
out having any preferred orientation.r| 

At the outer radius r = ro, the disc is expected to be 
truncated by a localized tidal torque associated with the 
m = 2 Lindblad resonance (Papaloizou & Pringle 1977; 
Paczyiiski 1977). The radius of the Lindblad resonance is 
unchanged by the tilt of the disc, and Larwood et al. (1996) 
have found that tidal truncation is only marginally affected. 
We apply the boundary conditions I x G — and v = 0. 
These imply that dl/dr = 0, and are appropriate if the 
truncating tidal torque is parallel to I, and if only angular 
momentum, rather than mass, is removed through the outer 
boundary. Again, these boundary conditions do not specify 
any preferred orientation for the outer edge of the disc. 

These boundary conditions are intended to be close to 
those used by Wijers & Pringle (1999), which were, however, 
implemented in an explicitly grid-dependent manner. 

2.3 Non-dimensionalization 

At least three reasonable choices of units are possible when 
implementing the equations numerically. First, the equa- 
tions could be implemented directly in CGS units. This 
might be appropriate when studying a particular system 
with well determined parameters, and has the advantage 
that the results can be compared immediately with the ob- 
served properties. In practice, however, many critical pa- 
rameters are poorly known, even for the intensely studied 
system Her X-1. Moreover, this direct approach obscures 
the scaling relations that exist between different systems. 

■f If the disc is terminated by a magnetosphere, a preferred orien- 
tation may well be introduced, but we neglect this possibility. In 
fact, recent analyses of the pulse shape changes in Her X-1 (Blum 
& Kraus 2000; Scott, Leahy & Wilson 2000) do not indicate that 
the inner disc is forced into alignment with the spin axis of the 
neutron star, as might have been expected. 
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Alternatively, a system of units could be adopted that 
is based on the geometry of the binary orbit. The radius 
and angular frequency of the binary orbit would be taken 
as basic units. This approach would be appropriate when 
studying a system in which the properties of the binary are 
well constrained, but perhaps the mass accretion rate and 
luminosity are less well known. 

In this paper we adopt a third system of units which is 
based on the accretion rate and the properties of the cen- 
tral object. This is particularly suitable for studying the 
radiation-driven instability in isolation, and the effect of the 
tidal torque can be added as an external perturbation. This 
will be important, for example, in order to demonstrate that 
the radiation-driven instability can naturally explain retro- 
grade precession even without the help of the tidal torque. 
This choice of units also has the major advantage of remov- 
ing the dependence on the mass accretion rate, one of the 
most poorly known parameters, in the dimensionless equa- 
tions. 

We therefore express r (and ri, To, re and rb) in 
units of GMi/c? (half the Schwarzschild radius of the 
central object), Q, in units of (? jGMx^ G in units of 
{M/2-k) {GMi /c) , J in units of (M /2n) (GMi /c^ ) , E in units 
of C~'^^'^°{M/2-Ky^'^°{c^/GMiy/^ and the modal or pre- 
cession frequencies {ui or Up, introduced below) in units of 
C^/i''(Af/2^)3/iO(GMi)-3/2c5/^ 

This simplifies matters such that the angular velocity is 



n 



-3/2 



the factor Cx drops out of equation (M), and the 



radiation torque becomes (cf. equation A15) 



6r 



■I X a. 



(15) 



where the dimensionless parameter e = Lt/Mc? is the ac- 
cretion efficiency. The tidal torque becomes (cf. equation 

b3) 



3Er^ 



Ttide = -/tido^^(e. ■ 0(ez X Z)(l + ■ ■ •), 
where 



/tide = qC~'l'\Ml2^r'l^\GM^cfl^ 

is a dimensionless parameter, with g = M2/M1. 



(16) 



(17) 



3 STABILITY OF A STEADY, FLAT DISC TO 
WARPING 

The first dynamical problem we consider is the linear sta- 
bility of an initially flat disc to warping. This problem was 
first considered by Pringle (1996) within a WKB approxima- 
tion, and later by Maloney, Begelman & Pringle (1996), Mal- 
oney & Begelman (1997) and Maloney, Begelman & Nowak 
(1998), who examined the general properties of solutions of 
the linearized equations and only later considered the effect 
of an outer boundary condition. Our investigation differs in 
that we consider a disc of given size, compute the discrete 
complex frequency eigenvalues of its bending modes, and 
thus determine whether it is stable or unstable. Moreover, 
our treatment of the hydrodynamics is based on the new 
evolutionary equations of Ogilvie (2000) which are derived 
directly from the basic fiuid-dynamical equations under the 
simple hypothesis of an isotropic alpha viscosity. Finally, we 



also consider the effect of self-shadowing on the linear sta- 
bility problem. In some previous work this was described 
as an intrinsically non-linear effect, but here we are con- 
cerned with perturbations involving tilt angles that exceed 
the small variations in H /r across the disc. 



3.1 Solution for a steady, flat disc 

For a steady, flat disc in the binary plane we have 

l = e,, G = G,e,, (18) 



^=0 ^=0 

dt ' dt 



and Trad = Ttidc = 0. The mass equation (Gl) can be inte- 
grated to give 

r«E = -|^H(re-r), (19) 

2n 

where H is the Heaviside function, and we have noted that 
V = at r = To- The radial velocity can then be eliminated 
and we are left with the angular momentum equation (H) in 
the form 



M d/i ,, , , 1 dGz 

Hire — r) = . 

2nr dr r dr 

With Gz = at r — ri, the solution is 

M 
Gz = --—(h-hi), n < r < re, 
Ztt 

M 



(20) 



re < r < ro. 



(21) 



where hi = h{ri), etc. 

From the general expression (|7|) for G, we have for a 
flat disc 



G. = QioXr'n' 
and 



(22) 



(23) 



Thus, given the parameters Mi, ri, Vc, ro, M , a and Cx, the 
quantities Gz,1, E and v can be obtained explicitly. 



3.2 Linear perturbations 

Now consider linear Eulerian perturbations of the basic 
state, denoted by a prime. Then 



— — + -—-(rv 1^ + rvl^ ) = 
at r or 

and 

— - + -Trirv'L + rvL') = -— - +T', 
at r or r or 

where 

L' = T.'hez + T,hl'. 



(24) 



(25) 



(26) 



It is sufficient to consider the x- and y-components of equa- 
tion (Eq), which read 

Ot r or r or 



ot r or r or 



(27) 
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Also 

rlf' 



dr 



or 



3o2< . ^ ^ 3^2 9^^ 



G'y = QwTrn^l'y + QaoXr'n"— ^ + QgoXr-'t^"— ^(28) 

Introduce a complex notation (Hatchett, Begelman & 
Sarazin 1981) in which 



W^C+il'y, G = G'^+iG'y, 



Then 

^,dW IdG ^ Mh^,^^, 

M /, dW d/i„A ,,, 



T = T^+ir^. 



(29) 



with 



G = QioXr^t^^VK + Oio^r^n^^ 



9VF 



(30) 
(31) 

2a(2 + ia) ' ^^^^ 

In this linear problem, only the first terms in the expansions 
of the Qi (equation UW are required. There is therefore no 
distinction at this stage between the models of linear and 
non- linear hydrodynamics defined in Section 2.2.1. 

If self-shadowing is neglected, the linearized radiation 
torque is (from Appendix A) 



where 

Qw = Q20 + iQao 



1 + 2ia + 6a^ 



U 



.dW 



.12ncJ dr 
Including self-shadowing, we obtain instead 



(33) 



-^ rad 



(Y^)i{[5l(emin + ^,0)-pi(e,nax,0)] 

dW 



+i [ff2(6lmin + VT, 0) - g2(6'max, 0)]} 

with 

51 (61,0) = i(6i + cos sine), 

TV 

92{0,O) = isin^e. 



dr 



(34) 



(35) 



The shadow angle 9s given by equation (A25) can be ex- 
pressed in terms of W as 



1" /r,r rrr-. dW 

= 2 + ^^'S(^ -W)- arg -^, 



(36) 



where W is the tilt variable of the ring that shadows the 
ring under consideration. 

The linearized tidal torque is (from Appendix B) 

3GM2T.r^ 



Ttidc = 



4r3 
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v2 175 
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\rh J 


' 64 
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X 



This can also be written as 
GMi'Sr 



iW. 



Ttid 



< 



'^/2 [rj 



iW, 



(37) 



(38) 



f,(m) 



where b\ is the Laplace coefficient of celestial mechanics. 
We may then seek a normal mode in which 



dW 
dt 



icoW, 



(39) 



where lo is the complex frequency eigenvalue. The problem 
then reduces to a complex linear eigenvalue problem involv- 
ing a second-order system of ODEs. It is convenient to use 
W and G as the dependent variables. 

The inner radius is a singular point since 1=0 there. 
Here we may apply the arbitrary normalization condition 
I¥ = 1. We also have dW/dr = 0, G = and dG/dr = 
— {M/2n)dh/dr for the regular solution. We then integrate 
from r = n to r = re. 

Owing to the presence of a (5-function in the angular 
momentum equation, G is discontinuous at r = re, while 
W is continuous (Maloney, Begelman & Nowak 1998). In a 
standard notation, the discontinuity in G is 



\G] = lim G 



Mh 

2n 



w. 



(40) 



It is trivial to implement this jump condition and integrate 
further from r = re to r = ro. 

At the outer radius we have the boundary condition 
dVy/dr = 0. Thus one complex quantity (co) must be guessed 
to start the integration, and one complex condition must be 
met at the outer boundary. Such a problem is readily solved 
numerically by Newton-Raphson iteration. 

3.3 Numerical investigation 

In his initial investigation of the radiation-driven instabil- 
ity, Pringle (1996) gave a simple approximate criterion for 
instability of a disc to warping. The disc should be unstable 
when 

^ ^ ^, (41) 

rs e^ 

where rs is the Schwarzschild radius of the central object 
and rj = V2/i^i is the ratio of effective viscosities. Under the 
assumptions of this paper. 



»7 = 



3Q20 
Qio 



2(1 + lo? 
a2(4 + cP- 



(42) 



in linear theory (see Papaloizou & Pringle 1983; Ogilvie 
1999). The approximate criterion can be rewritten for our 
purposes as 



^b 



> 



l^-K^rf' rb 



GMxjc^^ £2 ^^- (43) 

The critical radius depends strongly on e and especially on 
a. Interestingly, it depends on the luminosity of the central 
object Li, only through the accretion efficiency e and not 
through the mass accretion rate. Based on the properties of 
a low-mass X-ray binary with g ~ 1 (e.g. Her X-1), we take 
»-o/rb = 0.3 (cf. Papaloizou & Pringle 1977; Paczytiski 1977) 
and re/rb = 0.09 (cf. Lubow & Shu 1975; Flannery 1975). 
We also take n = QGMi/c . The consequences of varying 
q and/or of assuming the mass input to occur at the outer 
radius are examined in Section 5.1 below. 

These considerations lead us to choose rb as the basic 
control parameter. It is also the quantity that varies most 
significantly and obviously between different systems. Two 
other important parameters that affect stability are a and 
e, which are not well constrained, but it is difficult to argue 
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why they should vary substantially between different sys- 
tems. Our approach, therefore, is to determine numerically 
the minimum binary radius for instability as a function of 
a and e and to compare this with equation (|43|). We then 
examine the effect, if any, of a tidal torque. 

In the absence of the radiation and tidal torques, the 
disc possesses a discrete set of bending modes with 0, 1, 2, . . . 
nodes in the eigenfunction 'Re{W) (say). These modes are 
all damped [Im(Lj) > 0] by viscosity and typically precess 
retrogradely [Re(w) < 0] as a result of the internal torque 
component with coefficient Q3. Mode is exceptional in that 
it is almost neutral {u} « 0) in comparison. This mode would 
be the trivial rigid-tilt mode if the system possessed perfect 
spherical symmetry. However, because the mass input oc- 
curs in a definite plane, the full rotational symmetry of the 
system is broken and mode becomes non-trivial. 

When the radiation torque is gradually introduced (e.g. 
by increasing e continuously from 0) the frequency eigenval- 
ues u} move and the number of nodes in the eigenfunctions 
may change. We do not attempt a proper classification of 
the modes in the general case, but continue to label them 
by continuity with the simple limit of zero external torque. 

In principle, to determine stability, one must check the 
imaginary parts of all the frequency eigenvalues of the disc. 
In practice, the higher-order modes can be neglected because 
they are strongly damped by viscosity. The eigenvalues are 
usually ordered such that mode becomes unstable first as 
rb is increased, and the precession is retrograde. Under some 
circumstances, as described below and also in Section 5.1, 
mode 1 can become unstable first, and the precession is then 
prograde. 

It is also possible to simplify the procedure by solving 
directly for a disc possessing a marginally stable mode. This 
is done by restricting ui to be real but treating Tb as a second 
eigenvalue. One then has two real eigenvalues instead of one 
complex eigenvalue, and Newton-Raphson iteration can be 
applied in two real dimensions. 

The results are shown in Figs 1 and 2. Fig. 1 shows the 
variation of the minimum binary radius for instability with 
a and e, in the absence of a tidal torque. It can be seen 
that the effect of self-shadowing is to increase the critical 
radius by about a factor of 2. The approximate criterion 
( [43I) reproduces quite accurately the steep scalings with a 
and e, but systematically overestimates the critical radius. 
The reason for this is probably that Pringle (1996) estimated 
the wavelength of the unstable mode as ;$ ro, whereas in 
fact the effective wavelength of the modified rigid-tilt mode 
is somewhat longer than this. More accurate estimates of 
the linear stability criterion were made by Wijers & Pringle 
(1999). 

Fig. 2 shows the effect of a tidal torque on the lin- 
ear stability. The torque acts against the instability, making 
the critical radius larger. Note that, with a significant tidal 
torque, mode 1 becomes unstable first when a or e is large, 
and this results in prograde precession. 



can be obtained only by a direct integration of the time- 
dependent equations. However, we can use the principles 
of bifurcation theory (e.g. looss & Joseph 1980) to try to 
anticipate and categorize this behaviour. 

When the control parameter r-b is increased through 
the critical value rcrit for linear instability, a branch of non- 
linear solutions bifurcates from the solution representing a 
fiat disc. Owing to the symmetries of the problem, these 
non-linear solutions are steady in a rotating frame of refer- 
ence, although the angular velocity of that frame remains to 
be determined as an eigenvalue. These solutions therefore 
represent steadily precessing discs, which are the simplest 
possible non-linear outcome of the instability. 

Close to the point of bifurcation, these solutions are 
very little warped and nearly obey the linearized equations. 
The shape of these solutions is predicted by the eigenfunc- 
tion of the marginally stable mode of linear theory, and the 
precession rate of the disc corresponds to the frequency of 
the marginally stable mode. 

The branch of non-linear solutions may extend initially 
either into rb > rcrit (a supercritical bifurcation) or into 
rb < rcrit (a subcritical bifurcation). It may then turn 
around, undergo further bifurcations to more complex solu- 
tions, or simply terminate (since the equations being solved 
are not globally analytic). This behaviour depends on the 
non-linear terms in the governing equations. 

Our approach is to solve directly for the non-linear so- 
lutions representing steadily precessing discs. Not only are 
these the simplest possible outcome of the instability, they 
are also appropriate for explaining those observed systems 
that display a regular precession. They can be found by 
solving a non-linear second-order ODE eigenvalue problem, 
much more quickly and accurately than by using a time- 
dependent method. Both stable and unstable solutions can 
be found, and the network of branches of solutions provides 
a framework for understanding the more complex non-linear 
behaviour that can result. 



4.1 Analysis 

For a steady disc precessing about the binary axis with an- 
gular frequency uip, we may set 



dT. ^ dL ^ 

-=0, — =a.pe.xi. 



(44) 



The problem is then reduced to solving a set of ODEs with 
LOp as an eigenvalue to be determined. 

The mass equation (0) may be integrated as before to 
give 



M 

ruS — —■T—H{rc — r) 

27r 



(45) 



4 STEADILY PRECESSING DISCS 

Having determined the conditions for instability of a fiat 
disc, we now consider the non-linear solutions that result. 
In general these solutions may be arbitrarily complex and 



The angular momentum equation (pi) then becomes 



IdG ^ Mh^, 
ujpB, X L = --—+T+ - — S{r - rc)(ez - I) 
r dr Ittt 

M d{hl) „, 



(46) 
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Figure 1. Critical binary radius for linear instability of a fiat disc to warping, in units of GMxji? . Left: variation with a at fixed e = 0.1. 
Right: variation with e at fixed a = 0.3. Dashed line: numerical calculation neglecting self-shadowing. Solid line: numerical calculation 
including self-shadowing. Dot-dashed line: approximate formula (fes). The tidal torque is not included. 
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Figure 2. Effect of a tidal torque on the critical binary radius. The lowest curve in each panel corresponds to the solid line in the 
corresponding panel in Fig. 1. The other curves show the effect of an increasing tidal torque, with (from bottom to top) /tide = 10* i 
2 X IC*, 3 X 10^. In some cases mode 1 (prograde) becomes unstable before mode 0, leading to a kink in the curve. The dashed segment 
is then thecontinuation of the marginal stability curve for mode 0. A typical value of /tide ~ 2.3 X IC* can be estimated, based on 
equation (jl^) with qK,\,CxK, 4.4 X lO^^ cm20/7g-3/Ts-i2/7_ f;^ ~. lO^^ gs"! and Mi ^ IAMq. 



It is most natural, when integrating the angular momen- 
tum equation, to use I and G as the dependent variablesPI 
In order to implement the equations numerically, one must 



S In fact the three components of I are not independent since it 
is a unit vector. However, the condition |J| = 1 may be used as 
a check on the accuracy of the numerical method. This condition 
was satisfied to within 10"^" for all solutions obtained. 



be able to determine both E and dl/dr from I and G. One 
may proceed by writing 



d' 1/1 
r-^ = mm, 



(47) 



where m is a unit vector satisfying Im = 0. Let n = Ixm, 
so that {l,m,n) is a right-handed orthonormal triad. Then 



G = Xr n^ (Oi l + Q2\i'\m + Q:i\ii\ n) 



(48) 
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Now construct the vector 
/ X G 



V 



I G 



— \i,\n-—Wm. 



Then 



f2+Ql 



\M 



(49) 



(50) 



is an equation to be solved for |?/)|. In the approximation of 
Unear hydrodynamics this is trivial. In the non-linear model 
it is not so, and one must be aware of the possibility of 
multiple roots, or of no root.pi 

Having obtained j-i/;|, one determines the coefficients Qi. 
Then 



1 = 



I G 



(3ir2n2' 

in conjunction with equation 



yields E. Now 






ttW'^^ 



which allows n to be eliminated, giving 



d[ 

dr 



Qi 



?i + ' 



(QiV +Q2I xv) 



(51) 



(52) 



(53) 



as required. 

At the inner radius of the disc, owing to the axisymme- 
try of the problem, we may set £y = there without loss of 
generality, and then I = (sin/3, 0, cos/3), where /3 is the in- 
clination of the inner disc. Therefore two real quantities, ujp 
and /3, must be guessed when starting the integration. We 
then integrate the angular momentum equation from r = ri 
to r — re- 

Both the surface density and the derivative of the tilt 
vector are discontinuous at r = re, where 






(54) 



while I itself is continuous. 

At the outer radius, we apply the boundary condition 
i X G = 0. This consists of two independent constraints on 
the solution (e.g. the x- and y- components of the equation 
I X G — 0, provided that the outer edge is not inclined 
at exactly 90°), which are sufflcient to determine the two 
unknown quantities uip and (3 by Newton-Raphson iteration. 

A similar non-linear eigenvalue problem was posed by 
Schandl & Meyer (1994). In their model, the external torque 
is due to a coronal wind. However, an earlier and incorrect 
form of the angular momentum equation was adopted, and 
the actual solution of the eigenvalue problem was not given. 



'I Calculations suggest that the right-hand side of equation (BOl) 
is a monotonically increasing function of |'0| under most circum- 
stances, and therefore its inverse is unique. An exception arises 
when a is small {a ^ 0.05), since then Qi changes from nega- 
tive to positive as |^| increases through the value |V'| = IV'Icrit ~ 
\/24a. In that case the inverse function is double-valued, having 
one value less than li/^lcrit and one value greater. One can imme- 
diately distinguish between these possibilities, since the condition 
X > implies that sgn(Qi) = sgn(I ■ G). However, this may not 
be relevant because the disc will typically become unstable to 
perturbations under these circumstances (Ogilvie 2000). 



4.2 Numerical investigation 

The easiest way to find non-linear solutions is to set Vh 
close to a value corresponding to marginal stability of any 
mode, and to search for a solution with a small value 
of P and with Up close to the frequency of the marginal 
mode. The branches of solutions can then be followed quasi- 
continuously as r^ is varied. 

The results of this calculation are shown in Fig. 3, where 
we have taken a = 0.3 and e = 0.1. The tidal torque is not 
included. We compare the results of a simplified model, in 
which self-shadowing is neglected and linear hydrodynamics 
assumed, with the full m.odel. The eigenvalue /3, which is the 
inclination of the inner disc, is plotted as an indication of 
the 'amplitude' of the solution. 

In the simplified model, a branch of solutions appears 
at precisely the point of marginal stability and rises very 
steeply out of a subcritical bifurcation. We label this branch 
because it results from the instability of mode in the 
flat disc. The branch reaches a maximum value of /3, turns 
around sharply and then extends over a wide range of rb. 
Further branches appear, as rb is increased, at points where 
modes 1,2,... become unstable. The branches do not inter- 
act with each other; it should be remembered that /3 is only 
a projection of an infinite-dimensional solution space. 

In the full model, the bifurcation diagram is qualita- 
tively similar except that the bifurcations are now supercrit- 
ical. The values of rb corresponding to nrarginal stability are 
greater, as noted in Section 3.3. 

In each case the solutions of branch precess retro- 
gradely, while those of branch 1 precess progradely. This 
shows that retrograde precession can be obtained naturally 
even in the absence of a tidal torque, as found by Wijers 
& Pringle (1999) but contrary to the conclusions of Mal- 
oney, Begelman & Nowak (1998). The precession results 
partly from the radiation torque and partly from the in- 
ternal torque component with coefficient Q3. As explained 
by Wijers & Pringle (1999), when the mass input occurs well 
inside the outer radius, the tilt usually increases outwards 
so that the radiation torque causes retrograde precession. 



4.3 Stability of the steadily precessing discs 

Whether the solutions representing steadily precessing discs 
are relevant to observed systems depends on their stabil- 
ity. Theory predicts that the first branch to bifurcate from 
the flat solution (here, branch 0) will be initially stable if 
the primary bifurcation is supercritical, but unstable if the 
bifurcation is subcritical. All subsequent branches are ex- 
pected to be initially unstable. 

In the subcritical case the first branch to bifurcate from 
the flat solution is expected to regain stability at the turning 
point (saddle-node bifurcation) where the minimum value of 
rb is achieved. At some point along the branch, however, a 
secondary bifurcation may occur at which the solution loses 
stability to a more complex solution. These expectations are 
shown schematically in the insets of Figs 4-5. 

We have attempted to verify these expectations by de- 
termining the stability of the steadily precessing discs. We 
begin by rewriting the basic equations in a frame of refer- 
ence that rotates at the precession rate of the disc, and then 
linearize about the steady solution. The linearized equa- 
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Figure 3. Bifurcation diagram of non-linear solutions representing steadily precessing discs. Tlie inclination of the inner edge is plotted as 
an indication of the 'amplitude' of the solution. Left: simplified model. Right: full model. The various curves represent distinct branches of 
solutions produced at successive bifurcations of the flat solution (corresponding to the horizontal axis). The solid line, branch 0, appears 
at the point of marginal stability, where mode of the flat disc becomes unstable. The dashed line, branch 1, extends to larger values of 
/3 (not shown). Here a = 0.3 and e = 0.1; the tidal torque is not included. 



tions admit normal-mode solutions proportional to exp(i(.<ji) 
and we seek to determine whether any solution exists with 
Im(a;) < 0, indicating the instability of the solution. This 
method is a generalization of the analysis of Section 3 and 
the solutions obtained there can be recovered with the more 
general method. 

Owing to the technical difficulty of this calculation we 
do not describe it in detail here. In short, the results of 
the calculation agree precisely with our expectations: in the 
simplified model where the transition is subcritical, branch 
is initially unstable, regains stability at the turning point 
and then loses stability at a secondary bifurcation; in the 
full model, the transition is supercritical and branch is 
initially stable before losing stability at a secondary bifur- 
cation. Branch 1 is unstable. 

We have performed this calculation under two difi'erent 
assumptions concerning the luminosity of the central source. 
In the first case the luminosity is constant and given by L* = 
eMc? , while in the second case the luminosity is variable and 
given by 



L* = e(-27rr-t;S) 



(55) 



This more realistic prescription follows the suggestion of Wi- 
jers & Pringle (1999) that the luminosity should reflect the 
instantaneous mass accretion rate at the inner radius, rather 
than the steady rate supplied to the disc. Note that the dis- 
tinction between constant and variable luminosity does not 
affect any of our other results. 

For the simplified model and in the case of constant lu- 
minosity, branch is stable in the interval 1.06 ^ rb/lO'' Ss 
2.45. With variable luminosity, it is stable for 1.06 ^ 
rb/10^ ^ 1.38 (Fig. 4). Therefore, in accordance with the 
expectations of Wijers & Pringle (1999), the variable lu- 



minosity decreases the stability of the solution. For the full 
model including self-shadowing and non-linear hydrodynam- 
ics, the ranges of stability are 1.82 ^ rb/10^ ^ 2.79 and 
1.82 ^ rb/10® ^ 2.63 (Fig. 5) for a constant and variable 
luminosity respectively. 

The interval of stability is small, but non-negligible, cor- 
responding to a variation of about 50% in the control pa- 
rameter, and this is consistent with the results of Wijers & 
Pringle (1999). If robust, this implies that strictly steadily 
precessing discs should be relatively uncommon among X- 
ray binaries, and this is consistent with observations. The 
more complex solutions produced at the secondary bifur- 
cation would be quasi-periodic, i.e. periodic in a rotating 
frame of reference. The accretion rate on to the compact 
object would vary periodically in such a solution. Further 
bifurcations probably lead to quasi-periodic or chaotic so- 
lutions which could be found only by integrating the time- 
dependent equations. 



4.4 Properties of the stable non-linear solutions 

The solutions on the stable part of branch (in the full 
model) have a consistent and physically reasonable form. 
As the control parameter r^ is increased from the point of 
marginal stability, the solution rapidly acquires a tilted and 
warped shape. At rb = 2 x 10® (see Fig. 6) the inner and 
outer parts of the disc have inclinations of 12° and 37° , re- 
spectively, to the binary plane. As rh is increased further, the 
outer tilt remains close to 40° while the inner tilt decreases 
somewhat, as seen in Fig. 5. 

For a given solution, the tilt angle /3 is an increa sing 
function of radius while the twist angle 7 (see equation Al) 
increases from at the inner edge to a maximum of not 
more than 50° at an intermediate radius before declining 
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Figure 4. Region of tlie bifurcation diagram where stable, steadily precessing solutions are found in the simplified model. These are 
shown with solid lines. (Variable luminosity is assumed.) Left: precession frequency in dimensionless units. Right: inclination of the inner 
edge. Inset; schematic bifurcation diagram showing the expected stability properties of the non-linear solutions in the subcritical case. 
Solid and dashed lines denote stable and unstable solutions, respectively. 
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Figure 5. Same as Fig. 4 but for the full model. Inset: schematic bifurcation diagram showing the expected stability properties for the 
supercritical case. 



somewhat towards the outer edge. The extent of the shadow 
increases outwards, covering about half of the outermost an- 
nulus. The amplitude \ip\ of the warp reaches a maximum, 
typically 0.4-0.6, at an intermediate radius. The values of 
|V>| attained fully justify our taking into account the non- 
linear expressions for the coefficients Qi and the radiation 
reduction factors gi and (72. The inner part of the disc is 
very flat, however, and therefore passive, justifying our ne- 
glect of the regimes dominated by Thomson opacity and/or 
radiation pressure. 



Further geometrical properties are discussed in Sec- 
tion 5.4 below. 



5 DISCUSSION 

5.1 Stability of X-ray binaries to radiation-driven 
warping 

There has been some discussion as to whether X-ray binaries 
could or could not be subject to radiation-driven warping. 
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Figure 6. Example of a steadily precessing disc on the stable part 
of branch in the full model. Here a = 0.3, e = 0.1, rb = 2 X lO** 
and no tidal torque is included. The precession is retrograde, 
with ojp = —5.6 X 10~^. Lines of constant r and (f (defined in 
Appendix A) are plotted. The binary system is viewed edge-on, 
showing the tilt and twist of the disc and the extent of the shad- 
owed regions (black). The inner and outer radii are inclined at 
12° and 37°, respectively, to the binary plane. 



We found that the approximate criterion (143) was close to 
the one found with a more complete treatment including 
boundary conditions and self-shadowing. This criterion de- 
pends very strongly on the viscosity parameter (rcrit tx a"* 
for small q) and on the accretion efficiency (rcrit oc e~^) as 
can be seen from Fig. 1. In principle, it is always possible to 
find a combination of a and e for which a given system is 
stable, slightly unstable, or highly unstable. 

Estimating a in accretion discs is a long sought goal 
and the best results so far have come from comparing the 
light-curves of dwarf novae and soft X-ray transients to the 
predictions of the disc instability model. According to these, 
a is in the range 0.1-0.3 when these systems are in outburst 
and unlikely to be larger (e.g. Smak 1999; Cannizzo 1998). 
This probably also holds for persistent systems. Unfortu- 
nately the value of a cannot yet be reliably determined from 
numerical simulations of magnetohydrodynamic turbulence 
in accretion discs. While some local simulations suggest val- 
ues of 0.01 or less, there are important dependences on the 
mean magnetic field strength and on numerical details such 
as resolution and boundary conditions (Brandenburg 1998) . 
More recent, global simulations of cylindrical discs or thicker 
tori suggests values of 0.1 or greater (Hawley 2000). As for 
e, not all X-ray binaries host maximally spun-up black holes 
for which a maximum efficiency of 0.30 is reached (Thorne 
1974; although see Gammie 1999). For most systems a = 0.3 
and e = 0.1 are therefore plausible, if somewhat optimistic. 



Figure 7. Stability of X-ray binaries to radiation-driven warping 
{a = 0.3 and e = 0.1). rb is the binary separation in units of 
GM\/c^ and q is the mass ratio M2/M1. The upper two lines 
correspond to the first two bending modes when ?'add=''c. For 
q > 0.2, mode is unstable first (solid line) and mode 1 is shown 
as a dashed line. For g < 0.2, the situation is reversed. The lower 
solid curve corresponds to mode 1 for ?'add=''o. The critical rb 
for mode is then much higher and is not plotted. Also shown 
are the observed properties of a sample of systems (see Table 1). 
High-mass X-ray binaries are shown with circles, soft X-ray tran- 
sients with diamonds, other low-mass X-ray binaries as squares. 
Black hole candidates appear as open symbols and neutron star 
primaries as solid symbols. Crosses represent systems for which a 
long time-scale has been reported. 



5.1.1 Dependence on the mass ratio q 

We now examine the dependence of the stability criterion on 
the mass input radius radd and on the mass ratio q — M2/M1 
with a — 0.3, e = 0.1 and no tidal torque. As in Sections 3 
and 4, we first assume that the mass input radius radd is the 
circularization radius re. The ratios rc/rb and ro/rb are in- 
terpolated as functions of q from the tables of Lubow & Shu 
(1975) and Papaloizou & Pringle (1977). The variation of 
the critical radius for the first two bending modes is shown 
by the upper two lines in Fig. 7. For q ^ 1, rc/ro is almost 
constant and the critical disc radius varies only slightly. The 
critical binary radius for both modes therefore increases with 
increasing q as rb/ro. For q ^ 1, rc/ro increases and the crit- 
ical radius for mode becomes larger than what the approx- 
imate formula (|43) predicts. For q ^ 0.2, the prograde mode 
1 becomes unstable first, presumably generating a branch of 
stable, progradely precessing solutions. 

The assumption that radd = t'c rnay not be realistic 
when studying the stability of an initially fiat disc. In this 
case, the incoming stream from the companion will impact 
the disc at the outer radius ro, although some of the mat- 
ter may overflow to be incorporated further inside (e.g. Ar- 
mitage & Livio 1998). When radd ~ fo, mode becomes 
unstable only for very large rb and mode 1 is the first to 
become unstable. Furthermore, the critical radius at which 
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mode 1 is unstable is smaller than when mass is added at 
the circularization radius (Fig. 7) . This calculation does not 
take into account the fact that the mass added at ro has 
too little angular momentum for a circular orbit at that ra- 
dius. This would produce an extra torque which might affect 
the stability properties. If a disc can indeed be unstable to 
radiation-driven warping when mass is added at the outer 
radius, but stable when it is added the circularization ra- 
dius, then such a system might display warping cycles: an 
initially flat disc with mass input at its outer edge becomes 
unstable and tilts; mass input then moves towards re where 
the disc then becomes stable to warping and resumes its 
initially flat shape. We refer to this region in Fig. 7 as the 
'indeterminate instability zone'. The opposite situation was 
discussed by Wijers & Pringle (1999) who argued that, be- 
cause of the higher surface density, a disc with mass input 
at the outer edge would be more stable. Our results sug- 
gest otherwise but a more accurate treatment of mass input 
would be needed to reach a definite conclusion. 



5.1.2 Comparison with observed systems 

We have gathered data on a sample of high mass X-ray bi- 
naries showing 'super-orbital' modulation and on a sample 
of low-mass X-ray binaries including soft X-ray transientsll] 
(Table 1). All of the systems listed by Wijers & Pringle 
(1999) can be found here. The variability of the high-mass 
X-ray binaries in the sample is not due to periodic accretion 
in an eccentric orbit. There is enough evidence that these are 
disc-fed and have circular orbits (e.g. Bildsten et al. 1997; 
Ilovaisky 1984; Petterson 1978). They are all (except LMC 
X-3, see Section 5.3.2) above the critical binary separation 
for our chosen parameter values, due to the combination of 
their long periods and high q. Most of the low-mass X-ray 
binaries are below or, depending on the assumption for the 
mass input radius, slightly above the critical binary sepa- 
ration because of their shorter orbital periods. These con- 
clusions cannot be affected by the errors on the values of q 
and M\ : rh is weakly dependent on q and M\ for a low-mass 
X-ray binary of known orbital period so that the values we 
used would have to dramatically overestimate M\ and un- 
derestimate q to make all the systems unstable. This is not 
likely. 

An accurate treatment of mass input will probably yield 
a critical binary separation in between the two extreme cases 
we have considered and so should not modify qualitatively 
the above conclusions. Choosing a higher value of a and/or 
e would rapidly make low-mass X-ray binaries more prone 
to radiation-driven warping but, as argued above, there is 
nothing favouring higher values than a — 0.3 and e — 0.1. In 
addition, a more realistic model would take into account the 
tidal torque: although its effects are stronger when q is high 
(equation |l7[). Fig. 2 shows that the tidal torque significantly 
increases the critical binary separation (this is also found in 
the calculations of Wijers & Pringle 1999). Furthermore, 
the comparison with steadily precessing systems disfavours 



1 1 Radiation-driven warping can be important for soft X-ray tran- 
sients only during their outbursts. They are then close to steady- 
state with M{r) Ki constant so that we feel justified in including 
them here. 



Table 1. Parameters for the sample of X-ray binaries displayed 
in Fig. 7. The high-mass X-ray binaries appear first, followed 
by low-mass X-ray binaries (Her X-1) and soft X-ray transients 
(GRO J1655-40). Where there was no available estimate in the 
literature, the value we assumed is shown in brackets. We took 
Ml = 1.4Mq for neutron star primaries. Periods are in days, 
q = Mi /Ml, Ml is in Mq and the binary separation r^ is in 
units of GMi/c^. 



System 



Ml Tb/lO'^ 



Cen X-3 


2.090 


140 


17.0 


1.4 


6.7 


4U1907-f09 


8.380 


42 


12.0 


1.4 


15.3 


SMC X-1 


3.890 


60 


11.0 


1.4 


8.9 


LMC X-4 


1.408 


30 


10.6 


1.4 


4.5 


Cyg X-1 


5.600 


142 


1.70 


10 


1.9 


LMC X-3 


1.706 


99 


0.50 


5.0 


1.1 


SS 433 


13.10 


164 


[1.0] 


[10] 


3.0 


Her X-1 


1.700 


35 


1.56 


1.4 


3.1 


Cyg X-2 


9.844 


78 


0.34 


1.4 


8.0 


Sco X-1 


0.788 


62 


[0.7] 


1.4 


1.6 


AC211 


0.713 


37 


[0.7] 


1.4 


1.5 


GX 339-4 


0.620 


240 


[0.7] 


[5] 


0.6 


4U1957-I-11 


0.390 


117 


[0.7] 


[5] 


0.4 


4U1916-05 


0.035 


199 


[0.1] 


1.4 


0.2 


4U1820-30 


0.008 


176 


[0.1] 


1.4 


0.1 


Cir X-1 


16.55 


- 


[0.7] 


1.4 


12.2 


1746-370 


0.238 


- 


[0.7] 


1.4 


0.7 


1636-536 


0.158 


- 


[0.6] 


1.4 


0.5 


4U1626-67 


0.029 


- 


[0.07] 


1.4 


0.2 


GRO J1655-40 


2.620 


_ 


0.33 


7.0 


1.1 


Aql X-1 


0.800 


- 


0.30 


1.4 


1.5 


Nova Oph 77 


0.521 


- 


0.15 


4.9 


0.5 


Nova Vel 93 


0.285 


- 


0.14 


1.4 


0.7 


Nova Mus 91 


0.433 


- 


0.13 


6.2 


0.3 


GRO J0422-I-32 


0.212 


- 


0.11 


3.6 


0.3 


GS2000-I-25 


0.344 


- 


0.08 


8.5 


0.2 


A0620-00 


0.323 


- 


0.07 


7.4 


0.2 


V404 Cyg 


6.460 


- 


0.07 


12.3 


1.3 


Cen X-4 


0.629 


- 


[0.4] 


1.4 


1.3 


XTE J2123-05 


0.250 


- 


[0.4] 


1.4 


0.7 



a smaller critical binary radius than shown here as we shall 
see in the next section. We conclude that radiation-driven 
warping may only be relevant to those low-mass X-ray bi- 
naries with long orbital periods ^ Id (i.e. mostly soft X-ray 
transients) . 



5.2 Steadily precessing systems 

In Section 4 we have investigated whether radiation-driven 
warping could lead to a stable steadily precessing disc. We 
found such solutions existed only in a narrow range of pa- 
rameters. This probably explains why so few systems (Her 
X-1, SS 433, LMC X-4) show distinct, repeatable cycles al- 
beit with some variations. There is ample evidence that the 
cycles in these systems are not due to a variation of the in- 
trinsic luminosities (e.g. Margon 1984; Woo, Clark & Levine 
1995; Heemskerk & van Paradijs 1989). 

Most other 'super-orbital' periods do not appear clearly 
on periodograms and could be argued to be time-scales, 
modulations or perhaps beating between different long-term 
periods. When steady precession is possible, it is usually ret- 
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rograde which is consistent with the observations of Her X- 
1, SS 433 and LMC X-4 (Gerend & Boynton 1976; Margon 
1984; Heemskerk & van Paradijs 1989). For q — 1, we found 
the maximum binary separation for which steady behaviour 
could exist was about 2.6 in units of lO'' GM\/c? . Her X-1 
and SS 433, which have g « 1, lie very close to this value 
considering the uncertainties: from table 1 we get r^ = 3.06 
for Her X-1 and rb = 2.97 for SS 433. A slight decrease of 
a or e would provide better agreement, although again we 
emphasize that the correct tidal torque for each system has 
not been included. It is also rewarding that, in Fig. 7, the 
next closest source to the instability curve (for radd = »"c) 
happens to be LMC X-4. All other systems are further away 
above the curve or in the 'indeterminate instability' region. 

Because of the choice of units and the uncertainties 
in the different parameters of the observed systems, it is 
difficult to evaluate the precession period predicted by the 
model. In the range of stable solutions for q = 1, the preces- 
sion frequency ujp varies between 3.2xl0~^ and 7.0x10"^ 
(Fig. 5) in units of Cy^'\M /2TTf/^'-\GNh)-'^'^c'-/^ (Sec- 
tion 2.3). Using Ci = 4.4x 10^^ in CGS units (Section 2.2.1), 
this predicts precession periods for Her X-1 between 22d and 
47d with M ^ 4x 10^^ g s~^ . This estimate for M comes from 
assuming e = 0.1 and using the X-ray luminosity quoted 
by Choi et al. (1994), corrected for the distance estimate 
of Reynolds et al. (1997). It also agrees with the analy- 
sis of Cheng, Vrtilek & Raymond (1995). For Eddington- 
hmited accretion on to a WMq black hole with e = 0.1 
(i.e. M « 1.4 X 10^^ gs~^), we obtain precession periods be- 
tween 140d and 310d, consistent with the measured period 
of 164d for SS 433. However, this system may be accreting 
at a super-Eddington rate with low radiative efficiency. If 
e is significantly below 0.1 then it is unlikely to be unsta- 
ble to radiation-driven warping. Despite the apparent agree- 
ment, we do not feel warranted to make further quantitative 
comparisons because: (i) the predicted precession period de- 
pends on the poorly known mass accretion rate; (ii) the tidal 
torque almost certainly contributes to the precession rate, 
and (iii) we have not fully investigated the dependence of 
the precession period on the parameters of the model {q, a, 
e and radd). 

Qualitatively, the dependence of the precession period 
on Ml will result in longer precession periods for black-hole 
candidates by about an order of magnitude (e. g;. S S 433 
against Her X-1 above and Corbet & Peele 1997)Q Also, 
the frequency of mode at marginal stability, where it gen- 
erates the branch of steadily precessing solutions, is only 
weakly dependent on q for q > I (3.2-3.6x10"'^ in dimen- 
sionless units) and decreases to about 10~^ when q = 0.01 
where, in any case, very few systems are expected to show 
steady behaviour (see previous section). Although a more 
detailed analysis would be required, the range of steady pre- 
cession frequencies will not be much difi'erent when q > 1. 
This is consistent with the fact that LMC X-4 has similar 
properties to Her X-1 (Porb, Pong, Mi) except for q « 10. 



** This strengthens the interpretation of the 106d period in the 
nuclear X-ray source of M33 as being due to disc precession 
around a 10 Mq black hole (Dubus et al. 1997). 



5.3 Aperiodic variable systems 

A study of the behaviour of unstable systems would require 
a time-dependent method and our aim here is to provide a 
basic framework for interpretation. We have demonstrated 
in the previous sections that unstable systems close to the 
stability limit will undergo steady precession. As a system 
moves along the branch of solutions away from the stable 
region, it will evolve to more complex solutions, probably 
showing quasi-periodic behaviour before evolving towards 
chaos (Section 4.3). The exploratory calculations of Wijers 
& Pringle (1999) have shown such a sequence of behaviour. 
Here, the position of each system in Fig. 7 can be used to 
infer its behaviour. 



5.3.1 Warp-driven variability 

The systems that seem beyond the stable range for steady 
precession are Cen X-3, SMC X-1, 4U1907+09, Cyg X-2 and 
Cir X-1. SMC X-1 shows clear quasi-periodic cycles in the 
RXTE ASM data between 50 and 60 days that are consis- 
tent with increased absorption by a warped disc (Wojdowski 
et al. 1998). In the RXTE ASM, Cen X-3 shows a succes- 
sion of ON-OFF transitions that might be reproduced by a 
variable warp (Iping & Petterson 1990; Priedhorsky & Ter- 
rell 1983). Tsunemi, Kitamoto & Tamura (1996) could find 
no correlations between the pulse period history and the 
luminosity of this disc-fed pulsar. They concluded that the 
observed luminosity probably did not reflect the mass accre- 
tion rate, which is consistent if the variability results from 
a radiation-driven warp. 

Although in the unstable region because of its extremely 
long Porb = 16.6d, Cir X-1 has no reported long-term period- 
icity. Johnston, Fender & Wu (1999) proposed the neutron 
star is in an eccentric orbit which would make the present 
study irrelevant. This may also be the case for 4U1907-I-09 
for which a non-zero eccentricity is reported (Bildsten et al. 
1997). 

For q ^ 0.2, mode 1 becomes unstable first, imply- 
ing that systems there would most likely show prograde be- 
haviour. Cyg X-2 is close to this region and, indeed, the 78d 
variability would be consistent with a warped disc precess- 
ing progradely (Wijers & Pringle 1999; Wijnands et al. 1996; 
see also the discussion below of Paul et al. 2000). 



5.3.2 Marginally unstable or stable systems showing 
variability 

Some systems lie in the 'indeterminate instability' zone, 
amongst which are LMC X-3, Sco X-1, AC211 and Cyg X-1. 
The listed Pong for these are more likely to be time-scales 
than real periodicities. Recently, Wilms et al. (2000) pre- 
sented a long-term X-ray and optical study of the black 
hole candidate LMC X-3. The long-term spectral changes 
are not associated with changes in the column absorption 
as expected for a warped disc. The variations are easier to 
explain if the source of hard photons (the corona) is varying 
due to e.g. changes in the accretion rate. In Cyg X-1, which 
is right on the critical radius in Fig. 7, there seems to be ev- 
idence for both a warp (Brocksopp et al. 1999) and changes 
in the accretion rate or geometry (Nowak et al. 1999b). 

A modulation of the accretion rate through the disc can 
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occur on time-scales of 10-100 days when a viscously unsta- 
ble disc is irradiated by a constant X-ray flux (Dubus et al., 
in preparation) . This would also explain the time-scales and 
observations of GX 339-4 and 4U1957-I-11 which he well be- 
low the critical binary radius and, possibly, the very special 
cases of the ultra-short Porb systems 4U1916-05 and 4U1820- 
30 (which does not even appear in Fig. 7) where disc preces- 
sion can be completely ruled out. In this case modulations in 
the disc accretion rate cause the observed variability while 
in the systems well above rcrit it is quasi-periodic or chaotic 
warping that causes the variability. Nowak & Wilms (1999a) 
and Hakala, Muhli & Dubus (1999) have favoured a warped 
disc for 4U1957+11 but there is actually little observational 
evidence to support this. 

The recent study of Paul, Kitamoto & Makino (2000) 
comparing the long-term variability of Cyg X-2 and LMC 
X-3 is particularly interesting in this context. They argue 
that in both cases the observed variability is actually a 
combination of different more or less stable periodic compo- 
nents. This is consistent with the idea that Cyg X-2 could 
be switching between different branches of warped solutions. 
Contrary to LMC X-3, there is evidence in Cyg X-2 that the 
variations in luminosity are not due to changes in M (Ku- 
ulkers, van der Klis & Vaughan 1996). If these were to be 
due to varying obscuration then one would predict a harden- 
ing of the spectrum at low intensities as the softer photons 
are absorbed (or because the hard photons come from a lit- 
tle obscured corona). Paul et al. indeed find that the (5-12 
keV)/(3-5 keV) hardness ratio is anti-correlated to the in- 
tensity. But there is no correlation for the (3-5 keV)/(1.5-3 
keV) ratio and they conclude this does not support obscura- 
tion. Further observations at softer energies may be required 
to settle this. 



5.3.3 Soft X-ray transients 

We have argued above that systems in the 'indeterminate in- 
stability' zone could show some hysteresis-type behaviour, 
switching between warped and flat discs. Such behaviour 
may enhance variations in the mass transfer rate through 
the disc. This varying low-amplitude warp would not nec- 
essarily show up strongly as varying absorption, except for 
high-inclination systems. The long- Porb soft X-ray transients 
which lie in the 'indeterminate instability' zone (V404 Cyg, 
Aql X-1, GRO J1655-40, Gen X-4) could very well host this 
type of behaviour. 

The soft X-ray transients can very roughly be divided 
into two groups: one, associated with short orbital peri- 
ods, showing the classic fast-rise exponential decay (FRED) 
light-curve (e.g. A0620-00, GS2000-h25, Nova Mus 91, GRO 
J0422+32) and the other, with longer orbital periods, show- 
ing a wider variety of light-curves with plateau phases or 
erratic decays (see Chen, Shrader & Livio 1997). The FRED- 
type light-curve can be accommodated within the disc insta- 
bility model, which provides the framework for understand- 
ing the outbursts of soft X-ray transients and dwarf novae, 
if X-ray irradiation and disc evaporation in quiescence are 
taken into account (Dubus et al., in preparation). 

Yet this model does not provide an explanation of the 
variety of light-curves observed for systems with long Porb • A 
particularly interesting case is GRO J1655-40 where the long 
plateau during outburst and the observed anti-correlation 
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Figure 8. The accretion disc for the solution in Fig. 6 as seen 
from the compact object. The azimuth angle is equivalent to the 
precession phase. An observer at a high enough inclination i sees 
the disc (grey areas) periodically crossing his line of sight to the 
central source. 



between the optical and X-rays were explained by the com- 
bined effects of increased mass transfer from the secondary 
and variable warping (Esin, Lasota & Hynes 2000; see also 
Kuulkers et al. 2000, who show the dipping behaviour of 
the source changed during the outburst). Our study hints 
that long-period soft X-ray transients may be subject to 
such variable warping and this could partly explain their 
unusual outburst light-curves. Because of the dependence of 
the critical radius on Mi, neutron star transients are more 
likely to be unstable to radiation-driven warping and they 
should more often show deviations from the prototypical 
FRED-type outburst light-curves. 



5.4 Warping and irradiation heating 

The stable, steadily precessing discs in the full model (Sec- 
tion 4), which includes self-shadowing, all show a similar 
structure as discussed in Section 4.4. The shape of the disc 
when r-h — 2 X 10^ is shown from the viewpoint of the com- 
pact object in Fig. 8. The tilt of the outer disc is quite large 
(~ 40°) but Shakura et al. (1998) deduced similar values 
for the disc tilt in Her X-1 from RXTE observations. For SS 
433, the jet is inclined from its rotation axis by about 20° 
(Hjellming & Johnston 1981), also comparable to the values 
we obtain for /3 at the inner edge (see Fig. 5). 

An observer at an inclination i close enough to 90° will 
see a succession of ON-OFF states as the central source dis- 
appears behind the disc for certain phases. Obscrvationally, 
the X-rays are not fully eclipsed during the OFF states, sug- 
gesting the emission region is extended, very much like ac- 
cretion disc corona sources (Shakura et al. 1998). Further 
studies would need to assess whether such partial obscura- 
tion is compatible with the present model. Interestingly, the 
disc in a system seen edge-on need not always obscure the 
central source as can be seen from Fig. 8; but we would need 
to include the effects of a non-zero disc thickness to verify 
this. 

Irradiation by the central source can heat the outer disc 
enough to quench the thermal-viscous instability. This is 
particularly important in low-mass X-ray binaries where the 
standard disc instability model without irradiation heating 
would predict many more unstable systems than actually 
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Figure 9. Another view of tlie stable, steadily precessing disc 
shown in Fig. 6. The grey-scale shows variations of the coefficient 
C. The regions intercepting a large fraction of the flux are dark 
while those that intercept little or none (shadowed) appear white. 
Re-emission from the two arms occurs on opposite sides of the 
disc. 



observed (van Paradijs 1996). Irradiation is most important 
in the outer regions but these are completely self-shadowed 
in a flat disc (Dubus et al. 1999). If, however, the disc is 
warped then the outer regions may see the central source. 
Following Appendix A, the amount of irradiation received 
by each element in the disc is 



fir 



F n 



Aiir^ 



(56) 



where we have defined C = e,. ■ n for the warped disc 
(Shakura & Sunyaev 1973; Dubus et al. 1999). 

We show in Fig. 9 the variation of C in the disc. Dark 
areas are able to intercept a large fraction of the flux from 
the central source while white areas intercept little, or none 
if they are shadowed. The strongly irradiated regions form a 
two-armed spiral, although it should be understood that the 
re-emitted radiation from the two arms would be emitted on 
opposite sides of the disc. The average value of C over an an- 
nulus, taking self-shadowing into account, can be integrated 
like the radiation torque and this is shown in Fig. 10. In most 
of the outer region C is quite important (~ 0.1) and even a 
large albedo (the fraction of the flux that is scattered) would 
not prevent significant heating of the outer disc. For a solu- 
tion very close to the stability limit and thus corresponding 
to a much smaller tilt (about 2°), we find C ~ 5 ■ 10"''. 
Dubus et al. (1999) found that values of C of the order of 
10"'^ would be sufficient to heat the disc significantly and 
stabilize it against the thermal-viscous instability. A warp 
can easily provide such values of C even with a significant 
albedo. Taking the variations of H/r into account would 
probably not change these conclusions: H/r would be ex- 
pected to increase with radius, albeit with some variations 
due to the warping (equation 43 of Ogilvie 2000), so that 
even a partially shadowed outer disc would see enough of 
the X-ray fiux. 

While X-ray irradiation appears important for all low- 
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Figure 10. Average value of C (see Fig. 9) as a function of r, 
taking self-shadowing into account. 



mass X-ray binaries (van Paradijs & McClintock 1994) , only 
a fraction would be expected to be unstable to radiation- 
driven warping (Section 5.1). If warping is to explain how 
the outer disc can intercept the X-ray flux in all X-ray bina- 
ries then the warp must be produced differently (e.g. winds, 
Schandl & Meyer 1994; magnetic flelds, although this may 
only work close to the magnetosphere, Lai 1999). The re- 
sulting disc will probably also be able to intercept a large 
fraction of the flux. Warping might have interesting con- 
sequences for Doppler tomography and eclipse mapping; if 
irradiation heating dominates, the emission from the disc is 
clearly asymmetric (Fig. 9; see also Still et al. 1997, where a 
search for such features is attempted for Her X-1); viscous 
heating would also result in asymmetric emission since it is 
signiflcantly larger where the disc bends most (equation 44 
of Ogilvie 2000). 



5.5 Uncertainties and neglected effects 

Our model rests on a number of assumptions and simpliflca- 
tions. Some of these may have affected our conclusions and 
we list them here as suggestions for future investigations. 

(i) Regarding the hydrodynamics of the warped disc, 
the principal uncertainty is the modelling of the turbulent 
stress. Although our results here favour a viscosity param- 
eter Q Jj 0.1, this depends on our assumption that the un- 
derlying effective viscous process is isotropic. This difficult 
issue can be addressed only through numerical simulations 
of the turbulence (Torkelsson et al. 2000). 

(ii) Although studying the dynamical effects of radiation, 
we have not taken into account the thermal effects of one- 
sided irradiation on the vertical structure of the disc. This 
has been considered by Ogilvie (2000) and could result in a 
fractional correction to equation (bI). 

(iii) In our treatment of irradiation and self-shadowing, 
we neglected the non-zero thickness of the disc. 
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(iv) Our treatment of mass input to the warped disc is 
admittedly oversimplified (Section 2.2.2). 

(v) The correct inner boundary condition for different 
central objects is debatable (Section 2.2.3). 



6 CONCLUSION 

We have studied the radiation-driven warping of accretion 
discs in X-ray binaries. The latest evolutionary equations 
were adopted, which extend the classical alpha theory to 
time- dependent thin discs with non-linear warps (Ogilvie 
2000). We have also developed accurate, analytical expres- 
sions for the tidal torque and the radiation torque, including 
self-shadowing, that can be easily implemented in numerical 
calculations. 

(i) Using the complete set of equations, we re-examined 
the stability of discs to radiation-driven warping. We found 
that the critical binary separation is within an order of mag- 
nitude of the approximate criterion given by Pringle (1996) 
if the ratio of effective viscosities ri is given by equation (|42). 

(ii) Only the low-mass X-ray binaries with the longest 
orbital periods (Porb ^ Id) are likely to be unstable to 
radiation-driven warping. This could explain the unusual 
outburst light-curves of the long-period soft X-ray tran- 
sients. The disc-fed high-mass X-ray binaries are more likely 
to be unstable. 

(iii) We solved directly for several branches of non-linear 
solutions representing steadily precessing warped discs. We 
studied the stability of these branches to further perturba- 
tions and found that typically only one branch is stable, and 
that only in a limited range of parameters. In this case, the 
precession is usually retrograde. 

(iv) Discs that are beyond this range of parameters will 
probably show quasi-periodic or chaotic behaviour. 

(v) For q ^ 0.2, the radiation-driven instability may pro- 
duce progradely precessing warped discs. 

(vi) There may exist a certain range of parameters (the 
'indeterminate instability zone') for which a disc might cycle 
between being warped and being flat. 

(vii) Our results are sensitive to assumptions concerning 
the effective viscosity of the disc. If the turbulence acts on in- 
ternal shearing motions similarly to an isotropic effective vis- 
cosity, then values of a exceeding 0.1, and preferably closer 
to 0.3, are required. 

Further studies should consider individual systems, in- 
cluding the correct tidal torque on a case-by-case basis, and 
should make a more detailed comparison with observations. 
For systems in which periodic behaviour is not indicated, a 
time-dependent method should be used to model their vari- 
ability. 

The present results confirm the conclusions of ex- 
ploratory calculations by Wijers & Pringle (1999). Our 
methods, which are based on solving only ordinary differ- 
ential equations, should be seen as complementary to their 
time-dependent numerical approach. High-precision solu- 
tions can be obtained rapidly using our approach, and traced 
throughout the parameter space, but the non-linear dynam- 
ics can only be followed up to the point at which the steadily 
precessing discs become unstable to further perturbations. 



We find that showing the binary separation as a func- 
tion of the mass ratio can be a powerful tool to discuss the 
behaviour of X-ray binaries with respect to radiation-driven 
warping (Fig. 7). For the systems selected, there is a very 
persuasive consistency between the model and the observa- 
tions: for instance. Her X-1, SS 433 and LMC X-4 are very 
well explained by stable steadily precessing discs. However, 
most long-term variabilities observed in low-mass X-ray bi- 
naries cannot be associated with a radiatively-driven warped 
disc and might be due to modulations of the accretion rate. 
If warps are found to be ubiquitous in low-mass X-ray bina- 
ries, then it is likely that other driving mechanisms are at 
work. 
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APPENDIX A: RADIATION TORQUE 

In this section we derive an analytical expression for the ra- 
diation torque density Trad- Our discussion is based largely 
on the analysis of Pringle (1996). 



Al Geometrical considerations 

Consider the disc to be composed of a sequence of concentric 
circular rings. The unit tilt vector l{r,t) of the rings varies 
continuously with radius r and time t. In the notation of 
Ogilvie (1999), 



I — (sin /? cos 7, sin /3 sin 7, cos /3), 



(Al) 



where /3(r, t) and 7(r, t) are the Euler angles of the tilt vector 
with respect to a fixed Cartesian coordinate system {x, y, z). 
Introduce the dimensionless complex variable 

97 



Then 

|V| = 



(dl 



i^'^'^l^+^i;^""^ 



dl 



(A2) 



(A3) 



dl I, 



1 91 I, 
rlx — = \^ 



(A4) 



(A5) 



dlnr 
is a measure of the amplitude of the warp. In particular, 

cos /3 cos 7 cos X — sin 7 sin x ' 

cos f3 sin 7 cos x + cos 7 sin x 

— sin /3 cos X 

— cos /3 cos 7 sin x — sin 7 cos x ' 

— cos f3 sin 7 sin X + cos 7 cos x 
sin /3 sin x 

The position vector of a point on the ring of radius r is 
r = rsr, where (equation 7 of Ogilvie 1999, with 9 = 7r/2) 

cos /3 cos 7 cos (j} — sin 7 sin < 

cos /3 sin 7 cos (j) + cos 7 sin ( 

— sin P cos (j) 

is the radial unit vector and 4> the azimuth on the ringFTj 
Note that I ■ Er = 0, since points on the ring lie in the 
plane that passes through the origin and is orthogonal to I. 
Now (r, (j)) are (in general) non- orthogonal coordinates on 
the surface of the disc. The element of surface area is 

9r , \ I dr 

d4>^ 



(A6) 



dS 



{» 



(A7) 



which simplifies to 

dS = [I + Itpl cos(</f> — x) Cr] r dr d(/!>. 
The unit normal to the surface is then 

si -1/2 

X)\ 



Fi 1 I ; |2 2 

n — \1 + \tp\ cos 



(A8) 



[1 + |V|cos(<j!>-x)er]. (A9) 



TT This definition of (f> differs by it/2 from that of Pringle (1996). 
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A2 Radiation 

Let the disc be illuminated by isotropic radiation from the 
central source, of luminosity L*. Then the energy flux is 



F ^ 



Anr^ 



(AlO) 



except in regions of shadow. The power incident on an illu- 
minated surface element dS is 



47rr^ 



dP= |F-dS| = ■j-^\ip\\cos{(j)-x)\rdrd<j>. (All) 



Assume that the element absorbs the incident radiation and 
re-radiates it immediately and uniformly from the same side. 
Since the incident radiation is parallel to the radius vector, 
the element receives no torque from the absorption process. 
However, it receives a reaction force 

dF = -|— (±n) (A12) 

S c 

from the re-radiation process, where (±n) is the normal 
pointing away from the disc on the side that received the 
radiation. The correct sign is determined from the condition 
Br ■ (in) < 0, which implies 



dF = 



L* 



I cos(0 — x) '^^dr ( 



(A13) 



Gyrr^c 
The corresponding torque is 

r xdF = -r-^q(l + q'^y^^'\er xl)rdrd(f>, (A14) 

DTrrc 

where q = l-i/)] cos(0 — x)- The radiation torque density 
Trad (r, t) is defined such that 27rTrad r dr is the torque act- 
ing on a ring of radius r and of infinitesimal radial extent 
dr. Therefore 



127rrc 



■I X a. 



where 



a ^ — q{l + q ) ' erd(f> 



(A15) 



(A16) 



is a dimensionless integral. 

Assume at this stage that there is no self-shadowing, 
so that the integral is taken from (p = to (p = 27r, or, 
equivalently, from ijf> = xto<^ = x + 27r. Since e^ is a linear 
combination of cos <j} and sin 0, the basic integrals required 
to compute a are 



1 



2tt 



(l + X COS a) ^' cos ada — fix), 



l + x cos q) cos a sin ada — 0, (A-17) 



'0 

where 



/ W - ,,2(1 + ,2)1/2 [(1 + ^")E{k) - K{k)] . (A18) 

Here 



r/2 



,2.2 N-1/2 



K{k) = / {l-ksmay-'da, 
Jo 

/•t/2 

E{k) = / {I - k"^ sin a)^^^ da, 
Jo 



(A19) 



are the complete elliptic integrals of the first and second 
kinds, respectively, of modulus k — x/{l + x^)^' ^. Thus 

1 '•'^ 



/^ . 2^-1/2 

q(l + q ) cos 



/o 
and so 

We obtain finally 



2N-1/2 



sin (f) d(t> 



i2.c^(i^i)'x|;- 



Icosx/dV'l), 
|sinx/(|V'l), (A20) 

(A21) 

(A22) 



The factor /(|'!/'|) reduces the effectiveness of the radiation 
torque when the disc is significantly warped. It is a mono- 
tonically decreasing function of its argument, with /(O) = 1 
and f{x) --^ 4/(7ra;) as a; ^ oo. It has the Taylor series 



f(x)^l~-x +—X ~-—x +0{x 



8 



64 



1024 



(A23) 



A3 Self-shadowing 

Consider the ring of radius r, with tilt vector I, and consider 
an arbitrary ring of smaller radius f < r, with tilt vector 
I ^ I. The smaller ring obscures the radiation in the plane 
orthogonal to I and casts shadows at two diametrically op- 
posite points on the larger ring. These points lie in the plane 
orthogonal to I and are therefore defined by r ■ Z = 0. This 
determines their azimuth (j>s according to 

^ ^j; cos /3 cos 7 -I- ^„ cos /3 sin 7 - ^2 sin /3 /»„.n 

tan<^s = ^-^ 7^- ■ -■ (A24) 

£x sin y — £y cos 7 

However, it is more convenient to measure the azimuth with 
respect to the angle Xi i-e. 6 = (j> — x- We then find the 
covariant expression 



tan( 



dr 



''■''^ IT 
or 



(A25) 



Note that Oa is not uniquely defined since any multiple of 
n may be added to it. Consider a single branch, therefore, 
and follow 6s continuously as f varies from the inner radius 
of the disc up to r. Let Smin and Smax denote the minimum 
and maximum values obtained. Then the shadow on the ring 
considered consists of the two segments ^min < 6 < 6'max and 

6'min-|-7r < 6 < 6'max-|-7r. If theSC overlap (i.e. if 6'max - 6'min > 

n), the ring is completely shadowed and the radiation torque 
is zero. Otherwise, we must reconsider the integral a in the 
form 

a = - I q{l + q^)-^^^erde 



e„,i„ + 27r 



q(l + q')-'^^erde 

■K 

q(l + q^)-^'''erde. 



(A26) 



Smai 



Now the basic integrals required are 

— / {l + x cos a)^ ' cos ada — gi{6,x) 
^ Jo 
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e 



— I (1 + X cos a) ' cos a sin a da — g2(6,x), (A27) 
where 

g2{9,x) = -^ Ul + xY^' - {I + x^ cos^ 9)''^] ■ {A28) 
Here 



F{9,k) = / (l-fc^sin^Q)"^/^dQ, 
'o 

,2.2 xl/2 



E{9,k) = / (l-ksinay'da, 



(A29) 



are the incomplete elliptic integrals of the first and second 
kinds, respectively, again of modulus k = a;/(l+a;'^)^". Thus 



a = [gi{9niin + TT,M) - gi{9ma^^ 
+ [g2(6',nin+7r, |^|)-g2(6'n 



dr 



dl 



rlx —. (A30) 
or 



Noting that gi{6 + tv,x) ~ gi{9,x) — f{x) and g2{9 + tt , x) — 
g2{9, x) = 0, we confirm that the previous result is recovered 

in the limit 6lmax — ^min -^ 0. 

The elliptic integrals can be evaluated quickly and ac- 
curately using Carlson's algorithms (see Press et al. 1992). 



A4 Scattering 

In reality, when X-rays from the central source are inci- 
dent on the disc, a significant fraction of the photons are 
scattered from the surface rather than absorbed. Wijers & 
Pringle (1999) suggested that this X-ray albedo might re- 
duce the effectiveness of the radiation torque. However, we 
consider that the torque due to scattered radiation will be 
very similar to that due to radiation that is absorbed and 
re-emitted, because in each case the photons leaving the disc 
have the same total energy flux as the incident radiation flux, 
and the distribution of momenta of those photons is peaked 
around the normal to the disc surface. The non-isotropy of 
the Thomson cross-section is quickly lost after a few scat- 
terings and the photons finally leaving the disc have paths 
almost perpendicular to the surface. The torque due to scat- 
tered radiation could only be much smaller if the net effect 
of the scattering process was to divert the momenta of the 
photons through a small angle, which would require specular 
reflection from a slightly warped surface. In conclusion, the 
albedo of the disc is probably irrelevant to radiation-driven 
warping. 



APPENDIX B: TIDAL TORQUE 

Assume that the companion star may be treated as a point 
mass M2 describing a circular orbit in the xy-plane. The 
position vector of the companion star is then 



r2 = rb(cosfibt,sinJlbt, 0), 



(Bl) 



The time-averaged gravitational potential due to the com- 
panion star is 



$2(r) 



GM2 

27r 



■ r2 



■d(J7bi). 



(B2) 



For \r\ < \r2\ we obtain the expansion 



$2 = - 



+ 



GM2 


(s2 _ 2z^) 
1+ ^ ,2 + 


5(5s*' - 


- 90s*z^ + 1208^/ 



3(38* - 24s2z2 + sz") 



lez" 



256rS 



64r4 



+ 



(B3) 



where s^ — x'^ + y^ . 

Consider a circular ring of radius r, centred on the origin 
and tilted at an angle /3 to the binary plane. Owing to the 
axial symmetry of $2 , we may take the tilt vector to be 



I = (sin/?, 0, cos/3) 



(B4) 



without loss of generality. The position vector of a point on 
the ring is then 



r = r(cos/3cos(j!), sin^, — sin /? cost 



(B5) 



where </> is the azimuth on the ring. It is clear from the 
symmetry of the situation that the net torque on the ring 
has only a y-component. Now the y-component of the torque 
per unit mass is 



d(^2 9<&2 

ax oz 

When averaged over <; 



(B6) 



{ty) = 



GM2 

_525rf_ 
'65536rg 



3r 



T2 sin 2/3 
8rf 



this gives 
45r'^ 



1024r 



-(2 sin 2/3-^7 sin 4/3) 



(5 sin 2/3 -f 12 sin 4/3 + 33 sin 6/3) + ■ 



.(B7) 



The tidal torque density Ttidc is obtained by multiplying 
by the surface density E.[2J Noting that Bz ■ I — cos/3 and 
Bz X I = sin/?ej,, we may write Ttidc in the covariant form 

3GM2T,r^ 



Ttid 



+ 



175 
512 



Art 

15 
32 



r-b 



{Bz ■1)(bz X I) 
[7(Bz-lf~3] 
[33{ez -Z)* -SO{ez 



rb 



■lf + 5]+- 



(B8) 



This series converges quite rapidly. The magnitude of 
the fourth term (not shown) in the series in braces is less 
than 0.0027 when r/rb < 0.3. This suggests that the series as 
truncated above is accurate to within a fraction of a percent 
for all the calculations in this paper. Of course, higher-order 
terms can easily be calculated if desired. 

If the torque is not averaged over the binary orbital 
period, it is found to contain an additional, oscillatory com- 
ponent with frequency 2flb- The oscillatory torque causes 
a small wobbling motion but is unlikely to be important 
otherwise, except in very thick discs, or discs much smaller 
than the standard tidal truncation radius (Lubow & Ogilvie 
2000). 



■fT See Ogilvie (1999) for a demonstration that the surface density 
of a warped disc is independent of (f>. 
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